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PREFACE 

Part  of  the  Project  RAND  research  program  consists 
of  basic  supporting  studies.  Hie  research  presented  In 
this  Memorandum  concerns  certain  mathematical  techniques 
as  applied  to  problems  of  radiative  transfer. 
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SUMMARY 

The  use  of  principles  of  invariance,  as  in  invariant 
imbedding  and  dynamic  programming,  leads  characteristically 
to  functional  equations  of  the  form 

f ( P )  “  Tn(fn(®(P))^'  n 

where  fQ(p)  is  known.  The  computational  solution 
proceeds  stagewlse,  with  f^  determined  from  a  knowledge 
of  fQ,  f2  detennined  by  f^,  and  so  on. 

In  general,  what  is  desired  is  the  transient 
behavior,  small  n,  and  the  steady— state,  or  asymptotic 
behavior  as  n  -*  oo .  In  a  number  of  significant 
processes— -radiative  transfer,  control  theory,  inventory 
theory,  and  Markovian  decision  processes  In  general— only 
the  asymptotic  results  are  of  interest.  Ihis  is  also  the 
case  in  the  application  of  gradient  techniques. 

In  this  paper,  we  shall  outline  the  application  of 
nonlinear  summabillty  techniques  to  radiative  transfer. 
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A  NOTE  ON  NONLINEAR  SUMMABILITY  TECHNIQUES 
IN  INVARIANT  IMBEDDING 


1.  INTRODUCTION 

The  use  of  principles  of  invariance,  as  in  invariant 
imbedding  111  and  dynamic  programming  [2],  [3],  leads 
characteristically  to  functional  equations  of  the  form 

(1.1)  fn+1(p)  “  Tn(fn(e(p)))*  n  -  0,1,2,..., 

where  fQ(p)  is  known.  The  computational  solution 
proceeds  stagewise,  with  f1  determined  from  a  knowledge 
of  f q,  fg  determined  by  f^,  and  so  on. 

In  general,  what  is  desired  is  the  transient 
behavior,  small  n,  and  the  steady— state,  or  asymptotic 
behavior  as  n  -*  00 .  In  a  number  of  significant 
processes — radiative  transfer  [  ll ,  control  theory  1 3] , 
inventory  theory  [43,  and  Markovian  decision  processes 
in  general  C  53 ,  [6]  — only  the  asymptotic  results  are  of 
Interest.  This  is  also  the  case  In  the  application  of 
gradient  techniques  [ 73 • 

In  some  cases  where  the  asymptotic  results  are  not 
of  primary  importance,  they  are  worth  obtaining  in  order 
to  test  the  accuracy  of  the  numerical  techniques,  since 
the  steady— state  solution  can  often  be  derived  by  other 
independent  means. 

It  is  evident  that  direct  step-by-step  calculation 
of  the  asymptotic  solution,  using  (1.1),  la  time  consuming. 
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Hence,  it  is  important  to  develop  extrapolation  techr- 
niques.  A  most  important  step  in  this  direction  is  the 
work  of  Shanks  [  8] ,  closely  related  to  the  QI>-algorithm 
of  Rutishauser  [9].  We  shall  discuss  below  some 
extensions  related  to  nonequally  spaced  observations  or 
calculations.  In  this  connection,  let  us  mention  the 
work  of  Kantorovich  and  Krylov  tio]  concerning  the 
improvement  of  convergence  of  Fourier  series  and  other 
types  of  sequences. 

In  this  paper,  we  shall  outline  the  application  of 
nonlinear  summabllity  techniques  to  radiative  transfer. 
In  a  separate  communication,  we  discuss  its  use  in 
connection  with  control  theory  [ ll] . 


2.  MOTIVATION  OF  METHOD 

The  fundamental  assumption  is  that  we  possess  an 
asymptotic  expansion  of  the  form 


(2.1) 


f  o.f 
n 


00 


+  2 
k 


ake 


-v 


as  n  -»  00 ,  or  that  we  can  sensibly  approximate  to  fn 
by  an  expression  of  the  type  appearing  on  the  right-hand 
side.  From  theoretical  considerations  in  Invariant 
imbedding,  dynamic  programming,  gradient  techniques  and 
elsewhere,  we  obtain  rigorous  demonstrations  of  this 
relation;  see  [12],  [13].  What  is  remarkable,  as  shown 
by  Shanks  [  8] ,  is  that  quite  accurate  results  are 
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obtained  even  when  an  asymptotic  relation  of  this  type 
does  not  hold,  and  even  more  remarkable  is  the  fact  that 
calculations  based  on  small  values  of  n  give  excellent 
estimates  for  f  . 

If  we  take  f  to  have  the  simple  form 
(2.2)  fn  °  foo  +  aebn, 
an  immediate  calculation  yields 


(2.3) 


n 

fn+l 

n+1 

fn+2 

(fn  +  f„*>  -  Sfn>' 


In  the  next  section  we  shall  describe  some  experiments 
with  this  simple  nonlinear  predictor.  As  mentioned 
above,  many  further  results  based  upon  more  sophisticated 
approximations  will  be  found  in  Shanks  [8]. 


3.  RADIATIVE  TRANSFER 

Let  parallel  rays  of  radiation  be  incident  upon  a 
plane-parallel  slab  of  finite  thickness  which  absorbs 

radiation  and  scatters  it  isotropically  (see  Fig.  1). 


Fig.  1 


Using  the  theory  of  invariant  imbedding,  we  obtain 
an  equation  for  the  diffuse  reflection  r(0,^,x)  which 
leads  to  a  feasible  computational  algorithm.  Numerical 
and  analytic  results  are  given  in  [  14]  for  a  representa¬ 
tive  set  of  input  angles  9,  output  angles  and 

thickness  x.  As  pointed  out  above,  it  is  of  Interest 
to  obtain  the  results  for  infinite  x  in  this  way  in 
order  to  compare  with  previous  results  of  Ambarzumian 
and  Chandrasekhar  obtained  in  another  fashion. 

If  the  albedo  for  single  scattering  is  0.9*  a 
thickness  of  6  mean  free  paths  is  required  to  saturate, 
i.e.,  to  obtain  a  reflection  coefficient  equivalent  to 
infinite  thickness  to  about  four  decimal  places.  For  the 
particular  case  of  input  and  output  angles  at  60°  to 
the  normal  we  calculated  the  reflection  coefficient  at 
thicknesses  of  0.00,  0.02,  0.04,  ...,  1.20  mean  free 
paths.  These  values  are  listed  in  the  following  Table  1, 
which  is  read  across  from  left  to  right  in  each  row  and 
from  the  top  row  to  the  bottom  row.  For  a  thickness  of 
1.2  mean  free  paths  the  calculated  value  is  O.23295887. 

Next  we  use  (2.3)  on  each  set  of  three  consecutive 
entries  in  Table  1  to  produce  the  entries  in  Table  2. 
These  are  predictions  of  the  limiting  value,  which  is 
0.272389. 

Using  the  predictions  of  Table  2,  we  can  use  the 
formula  once  again  to  produce  another  set  of  predicted 
values.  These  are  shown  in  Table  3. 
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We  see  that  we  may  predict  a  limiting  value  of  about 
0.2J,  which  is  quite  accurate  enough  for  many  purposes. 

Also  remarkable  is  the  fact  that  a  double  precision 
calculation  of  values  at  .005*  .010,  ...,  .300,  with 
r(ir/3,  ir/3 ,  *300)  =  .111595  (which  is  less  than  50  per 
cent  of  the  limiting  value),  predicts  a  limiting  value 
of  0.27. 

It  is  clear  what  a  great  saving  in  time  can  be 
obtained  in  this  way. 

Harriet  Kagiwada  carried  out  the  calculations  on  an 
IBM-7090.  Only  a  few  minutes  of  computing  time  were 
required . 

The  questions  of  which  predictor  formula  to  use, 
what  increment  in  thickness  to  employ,  and  how  many 
increments  to  use  in  order  to  predict  the  limiting 
values  most  efficiently,  are  still  open. 

4 .  TIME-DEPENDENT  PROCESSES 

As  indicated  in  a  previous  paper  [ 15] ,  time-dependent 
radiative— transfer  problems  may  be  resolved  computationally 
in  a  multistage  fashion.  First  we  use  invariant  imbedding 
to  obtain  a  set  of  nonlinear  partial  differential— integral 
equations.  A  Laplace  transform  then  reduces  these  to 
equations  identical  in  form  to  those  encountered  in  the 
time— independent  problem.  These  are  integrated  numerically 
for  appropriately  chosen  values  of  the  transform  variable. 
Finally,  we  use  a  numerical  inversion  of  the  Laplace 
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transform.  A  disadvantage  of  this  method  is  that  the 
values  of  the  desired  function,  say  u(t),  are  obtained 
at  irregularly  spaced  points  t^,  i  =  1,2,...,M,  and  it 
is  not  easy  to  make  t^  large.  Hence,  an  accurate 
extrapolation  method  would  be  quite  useful  in  this  case. 

In  place  of  asking  that  we  have  a  representation  of 
the  type  appearing  in  (2.1),  we  can  ask  that  u(t)  be 
approximated  to  by  a  function  w(t)  satisfying  the 
linear  differential  equation 

(4.1)  w^  +  b^w^1*”1)  +  •••  +  bnw  =  bn+1. 

The  unknown  constants  b^^  and  the  Initial  values 
w^(o)  =  c1,  i  =  0,1,..., n  -  1,  are  to  be  determined 
by  the  condition  that 

(4.2)  2  (u(t . )  -  w(t, ) )2 

i=l  1  1 

is  a  minimum.  The  numei’ical  solution  of  problems  of  this 
type  can  be  carried  out  quite  easily  using  the  techniques 
of  t 16] . 

5.  GRADIENT  TECHNIQUES 

The  general  idea  of  the  gradient  technique  is  to 
solve  an  equation  of  the  form  T(u)  =  0  by  imbedding  it 
within  the  solutions  of 


(5.1) 


du 

cftf 


T(u). 


The  solutions  of  the  original  equation  are  taken  as  the 
steady— states  of  (5*1) •  Since  in  this  case  only  the 
values  at  t  =  oo  are  desired,  nonlinear  summability 
results  will  save  a  good  deal  of  computing  time.  Combined 
with  quasilinearization  [ 17l ,  t 18] ,  very  accurate 
results  can  be  obtained.  Results  of  this  nature  will  be. 
presented  subsequently. 
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